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ABSTRACT 



We present a technique to treat systems with very many levels, like molecules, in 
non-LTE. This method is based on a superlevel formalism coupled with rate operator 
splitting. Superlevels consist of many individual levels that are assumed to be in LTE 
O ' relative to each other. The usage of superlevels reduces the dimensionality of the rate 

equations dramatically and, thereby, makes the problem computationally more easily 
treatable. Our superlevel formalism retains maximum accuracy by using direct opac- 
ity sampling (dOS) when calculating the radiative transitions and the opacities. We 



developed this method in order to treat molecules in cool dwarf model calculations 
in non-LTE. Cool dwarfs have low electron densities and a radiation field that is far 
from a black body radiation field, both properties may invalidate the conditions for 
the common LTE approximation. Therefore, the most important opacity sources, the 
molecules, need to be treated in non-LTE. As a case study we applied our method to 
carbon monoxide. We find that our method gives accurate results since the conditions 
for the superlevel method are very well met for molecules. Due to very high collisional 
cross sections with hydrogen, and the high densities of H2 the population of CO itself 
shows no significant deviation from LTE. 

Subject headings: molecular processes - methods: numerical - stars: atmospheres - 
stars: low-mass, brown dwarfs - radiative transfer - line: formation 



1. Introduction 

Frequently, the atmospheres of late type dwarfs are considered to meet the conditions in 
which the approximation of local thermodynamic equilibrium (LTE) can be applied. In LTE the 
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population densities of atoms or molecules can be described by the TE (thermodynamic equilibrium) 
Boltzmann-Saha distribution, leading to very simple calculations. LTE requires that the collisional 
rates are large enough to fully compensate for the deviations of the radiative rates from their 
Planckian thermodynamic equilibrium values. Any spectrum that is not a Planckian radiation 
field will drive a system into non-LTE unless the collisional rates are strong enough to prevent 
this. In the extremely cool atmospheres of M dwarfs, L dwarfs or T dwarfs, however, there are two 
important effects which can lead to deviations from LTE. The first is that the extremely low electron 
densities and consequently low collisional rates might not be sufficient to restore LTE. Collisions 
with other particles might be much less effective because of their much smaller relative velocities 
and their smaller cross sections. Investigations for stars like the sun showed that electron collisions 
with CO molecules are negligible compared to all other excitation and de-excitation processes (e.g. 
Thompson 1973). For cooler stars, the electron collisions are even less important for all molecules, 
not only CO. Thompson (1973) also investigated the effects of collisions with H, H2 and He on 
the CO molecule and found them to be important under certain circumstances. The very detailed 
study by Ayres & Wiedemann (1989) collected the various existing collisional cross sections of CO 
with H, H2 and He. They found CO to be in LTE in solar like and cooler stars due to the very 
high quasi resonant cross section with H. 

Other molecules might be less affected by atomic collisions. In general, values for collisional 
cross sections are only known very roughly through the formulae of Drawin (1961), Van Regemorter 
(1962) or Allen (1973). Without good data on collisional excitation cross sections we cannot be 
certain that they will dominate the population distribution in molecules. 

In addition, with decreasing effective temperature, the maximum of the energy distribution 
stays at roughly 1-1.2 ^m (Allard & Hauschildt 1995) since the strongest opacity sources, TiO in 
the optical and H2O in the infrared, leave only a small window between 1—1.2 /im with relatively 
little absorption. This means that not only the shape of the spectra are far from spectra of black 
bodies for T = T c g, but also that the maximum deviates strongly from the maximum of a black 
body and the radiative transitions "see" a much hotter local radiative temperature compared to the 
kinetic temperature. This circumstance is particularly important for cool stars which have broad 
molecular absorption bands compared to narrow atomic lines for hotter stars. Of course, for both 
cool and hot stars, the radiation field produces non-Boltzmann populations if it is diluted due to 
optically thin layers, regardless of its shape. 

With these different temperatures the radiative rates and the collisional rates will try to pop- 
ulate the levels differently. Non-LTE effects for atomic species such as Ti have already been 
investigated by Hauschildt et al. (1997a, 1996a). That means that the general conditions for non- 
LTE effects (sufficiently low collisional rates and a sufficiently non-Planckian radiation field) are 
met in M dwarf atmospheres. Since molecules are the dominating opacity sources in cool objects, 
deviations from LTE can have a significant impact on the atmospheric structure and the spectra. 

These reasons lead us to develop a method to treat molecules in detailed non-LTE. We want 
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to emphasize that it is quite well possible that certain molecules will be in LTE under certain 
atmospheric conditions (as discussed above for CO when dominated by collisions with atoms and 
other molecules). However, the aim of this paper is to demonstrate how it is possible to calculate 
non-LTE effects in molecules in an accurate fashion and to lay the foundations for calculating 
molecules that are not in LTE in cool atmospheres (e.g. molecules that are insensitive to collisions 
or molecules that form their spectra in a part of the atmosphere that cannot sustain the conditions 
for LTE). 

Molecules have many more levels as well as many more lines compared to atoms. Several thou- 
sand levels and several hundreds of thousand lines are minimal numbers for very simple molecules. 
Complicated molecules (such as H2O) will have orders of magnitudes more levels and lines. The 
latest line data available for H2O (Partridge & Schwenke 1997) and TiO (Schwenke 1998) contain 
330 and 140 million lines based on about 50 and 20 million levels respectively. Since the system 
of rate equations (equation (4), below) is a system with a rank of at least the number of involved 
levels, it is obvious that a straight forward treatment of molecular non-LTE will exceed the limits 
of even the largest modern parallel super-computers. 

The huge number of lines and levels in molecular systems is merely a technical difference 
compared to atomic systems. There is, however, also an important physical difference. Atomic 
processes only include excitation and ionization processes (and their inverse processes). In addition, 
molecular processes include dissociation and formation of the molecule through different possible 
channels. That means that a realistic model has to include not only the molecule itself but also 
take into account the constituting atoms and the equation of state has to be aware of a possible 
source or sink of atoms. The logical extension is "Non Local Chemical Equilibrium (NLCE)". As 
a first step, we address in this paper only internal non-LTE and do not account for dissociation. 
Neglecting NLCE is a good approximation for very stable molecules which are not easily dissociated 
since the kinetic and radiative energies in M dwarf or cooler atmospheres are too small to do so. 

So far only very little work has been done on general molecular non-LTE mainly because 
of the difficulties mentioned above. Usually, only very few levels close to the ground state are 
considered. CO, however, has experienced a number of investigations. The limiting factor for CO 
is the uncertainty in the quality of the very high collisional rates. A recent and important example 
is the work by Ayres and Wiedemann (Wiedemann et al. 1994; Ayres & Wiedemann 1989) who 
calculated CO in non-LTE in order to try to explain the presence of strong CO bands in the solar 
infrared spectrum, which cannot be explained by the high kinetic energies of the gas in the solar 
atmosphere. They only used 10 vibrational levels (which produce the optically thick lines) and 
the rotational levels were simplified in a fashion similar to what will be presented here, i.e. they 
assumed some LTE behavior between certain levels. 

Another example is the work by Kutepov et al. (Kutepov et al. 1997, 1991), who calculate the 
solar radiation in the earth's atmosphere through molecular bands including CO. They consider 
only three vibrational levels and started with so called vibrational LTE (which is also conceptually 
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very similar to the work presented here) but relaxed this assumption in their latest work. 

In the next section we provide the necessary background for this work. We review the atmo- 
sphere code used and the physics involved. In section (3) we present the superlevel method on 
which this paper is based on. In section (4) we apply our method to carbon monoxide. We use CO 
as an example and demonstrate the numerical implementation. We conclude our paper in section 
(5). 



2. Background 

This work is based on the general stellar atmosphere code PHOENIX described in Hauschildt 
et al. (1996b); Hauschildt et al. (1997); Hauschildt et al. (1997b); Baron & Hauschildt (1998); 
Hauschildt et al. (1999). 

The solution of the wavelength dependent radiative transfer equation is performed by operator 
splitting as described in Hauschildt (1992) where the mean intensity J is iteratively calculated via 

J ncw =A*5 I1CW + (A-A*)5 i d . (1) 

where S is the source function, A the lambda operator expressing the formal solution and A* a 
suitably chosen approximate lambda operator. 

One major feature is that phoenix calculates each wavelength point separately, i.e. phoenix 
steps through the wavelength grid and solves the radiative transfer equation at each wavelength 
point by only using data important for that wavelength point. This dramatically reduces the 
memory requirements for the code. Only the quantities that need to be integrated over wavelength 
are kept in memory. 

Another important feature to be mentioned here is the use of line lists which contain several 
hundred million spectral lines, phoenix has been especially optimized to handle this task by 
using dynamical opacity sampling (dOS) (see also Schweitzer 1999; Hauschildt et al. 1994). dOS 
dynamically discards the negligible lines for a particular atmosphere configuration and does not 
require precomputed tables. With this feature phoenix can calculate spectra with high resolution 
and a large number of spectral lines efficiently The superlevel method presented below will take 
great advantage of dOS as it will improve the accuracy and reduce the necessary approximations. 

In order to calculate the population densities of individual levels in molecules (or atoms) 
without the assumption of thermodynamic equilibrium, we need to solve the rate equations which 
balance all population and de-population processes for every level. 

The radiative absorption rates for individual levels due to absorption of photons are described 
by the Einstein Bi u coefficient 2 and the profile averaged mean intensity J u \ or by the cross section 



2 The indices used throughout the paper have the following meanings : I always stands for a lower level, u always 
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af^ and the mean intensity J : 

a roc 

Riu = B lu J ul = — ^ «r« S ^(A)AdA (2) 

The emission rates are the sum of induced emission and spontaneous emission. However, to simplify 
the rate equations we shall follow Mihalas (1978) and use a radiative emission rate coefficient R u i 
between individual levels 

Rm = y c [ <i (^T + J a(A)) e-^AdA (3) 

such that we can write the emission rate as n u (nj /n^)R u [. 

With these radiative rate coefficients and with the collisional rate coefficients C u \ = (ra*/n*)C/ u 
the rate equation for one particular level i can be written as 

J2 MRh + Cu) + n u ^(Rui + c iu ) = m (j2(Riu + c iu ) + Y ^(Ru + c H )) ■ ( 4 ) 

l<i u>i u u>i l<i 1 

The system of all rate equations is closed by the conservation equation for the particles and the 
charge conservation equation. For a multi level system, equation (4) is a non-linear system of 
equations with respect to the population densities since the radiative rate coefficients themselves 
depend on the population densities via the mean intensities. It is also non-linear with respect to 
the electron density via the collisional rate coefficients and the charge conservation equation. In 
the case of molecules that can dissociate and form through various channels the situation will get 
more complicated since system (4) will include not only the molecule in question but also all atoms 
and all possible molecules which have to be calculated simultaneously. Furthermore, the condition 
of particle conservation has to be exchanged for the condition of nuclei conservation. 

We use the same formalism as Hauschildt (1993) and solve equation (4) with rate operators 
and operator splitting techniques. Then the rates can be written by means of a rate operator [Rij] 
and a population density operator [rij] : 

R ij = [R ij }[n i ] (5) 

The radiative rates are functions of the mean intensity J which is usually expressed by the lambda 
operator and the source function J = AS. Therefore, the rate operator is also a function of the 
lambda operator. However, since the source function is non-linear with respect to the population 
densities, the lambda operator is substituted by the \& operator introduced by Rybicki & Hummer 
(1991) 



A(A) = (A) 



_(«A + cr\) 



(6) 



stands for an upper level and i or j are used where a distinction between upper and lower level is not possible or not 
necessary. 
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This will eliminate the full source function and introduce the emissivity rji u {\) instead, which can 
be expressed by the operator [-E(A)] by 

[E(X)][n]=J2mu(X) + fj(X) (7) 

Ku 

where 77(A) is the background emissivity. 

With these definitions, the absorption rates in eq. (2) can be expressed as 



47T 

[Rlu][n] = Vc 



00 

a^(X)E(X)XdX 





[n] (8) 



and the emission rate from eq. (3) becomes 

poo /OA„2 



[Rul][n] = h~cj a ° J (-JF + *(A)[^(A)][n]J e~ hc / XkT XdX (9) 

In order to solve the system of rate equations with an operator splitting iteration scheme, the 
rate operators are split in analogy to the approximate lambda operator (1) : 

Rij = [i%][n new ] + ([Rij] - [i?y)[n old ] = [i%][n new ] + [A%][n old ]. (10) 

The rate equations (4) can then be written in the form 

^ ] n l,ncw [Rli] [^new] + ^ ] re «,ncw ~ [R-ui] [ n ncw] 
Ki u>i u 

([AR ui ][n old ] + C iu ) 

Ki u>i u 

u>i Ki 

(^([AR.JM + C iu ) + ^([ A ^z]Kid] + C K )). (11) 



77/ ■ 

This equation can be solved iteratively and represents an operator splitting iteration analogous to 
solving the radiative transfer equation. 

In phoenix, equation (11) is solved after having solved the radiative transfer equation for all 
wavelength points. During the wavelength loop all the rates and operators are integrated using the 
newly calculated mean intensity. Convergence of the occupation numbers is obtained by iterating 
this calculation with updated occupation numbers (and a possible temperature correction applied) 
until overall convergence. 
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3. The superlevel formalism 

3.1. General concept 

For complex molecules, the dimension of eq. (11) becomes very large (see sec. 1). However, 
the memory requirements for a computer program that has to solve eq. (11) with a direct approach 
will scale with the number of levels squared. Therefore, depending on the molecular system and 
the available computer resources the problem of solving the rate equations becomes untreatable 
and needs to be simplified. One possibility to simplify non-LTE problems is to treat the system 
"partially" in LTE. This was done, e.g., by Kutepov et al. and Ayres and Wiedemann as mentioned 
above. They assumed LTE population within one vibrational state. 

A much more general approach was presented by Anderson (1989) who faced the problem 
of complex atoms and ions with a lot of levels (e.g. Fe II, Ni II etc.). He also grouped several 
levels into one common level and assumed LTE within such a group. However, the exact grouping 
was relaxed to "some similarity" of all levels within one group. Also, he derived methods how to 
calculate the opacity by such levels by means of opacity distribution functions. This idea was also 
the basis of the work by Dreizler & Werner (1993) and by Hubeny & Lanz (1995). 

With the rapid development of large computers it is no longer necessary to use superlevels for 
atomic systems. Hauschildt et al. (1997) and Hauschildt et al. (1996b), e.g., calculate large atomic 
systems directly in non-LTE. For molecules, however, the number of levels is orders of magnitudes 
larger and approximations are still necessary. Furthermore, these approximations are much better 
for molecules than they are for atoms. The typical energy differences between individual molecular 
levels are much smaller than the energy differences between atomic levels. This provides very good 
internal coupling and thermalization within one superlevel. 

We will only adopt the basic idea of superlevels, that is, we group many individual levels 
into one superlevel and assume LTE populations within one superlevel. The populations of the 
superlevels will then be calculated via the rate equations. However, the absorption coefficients 
and the transition rates between the superlevels will not be approximated. We will account for all 
individual lines by dynamically sampling the opacity and the transition rates. This distinguishes our 
approach from the original approaches mentioned above. There, it had been necessary to develop 
opacity distribution functions for the transitions between the superlevels. This is not needed here. 
As explained above, phoenix uses dOS and, therefore, a dynamical opacity sampling of transitions 
between superlevels is already done by design. The only purpose of the superlevel formalism is to 
keep the system of rate equations small. 

A superlevel 3 / is constructed out of a number of actual levels i and its population number 



3 In the following, the same nomenclature conventions apply as in section (2) with the addition, that capital letters 
like I, L or U are used for superlevels, and that lower case letters like i, u or I are used for actual levels. 
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density is defined by 

m = y^m. (12) 

iei 

For the population densities of the actual levels within one superlevel LTE is assumed and the 
Boltzmann equation can be used for a level iei: 

n % = nj^e-^ kT (13) 

where 

Z I = Y j me- E * /kT (14) 
iei 

is the finite partition sum over the superlevel considered. 

Eq. (13) can be used for both the LTE and the non-LTE quantities. With this we define the 
departure coefficients as in Menzel & Cillie (1937) and Mihalas (1978) as 

bi = — = — = bi, (15) 

n* nj 

where the n* are the occupation numbers calculated via the Saha Boltzmann equation and the 
actual non-LTE population density of the continuum. Also, in LTE equation (13) holds true for 
any level i £ I. Summing equation (13) over all levels i € J we get the ratio nj/nj of two superlevels 
to : 

4 = ^ (16) 

which is the analogous relation to the Boltzmann relation n*/n* = gi/gj exp(-EijfkT) for normal 
levels and is in fact only a more general form of equation (13). 



3.2. The partition functions 

Equation (16) can be written as 

n 4 = ^ (17) 
n*j n*j Zi 

Summing over all superlevels and adding all remaining levels not included in any superlevel yields 

^ = ^ £/ bjZj + Y,i0 L e W (-E t /kT) (18) 
This leads to the definition of the non-LTE partition function which can be written for superlevels 
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where <5 CO rr is the sum over all levels not included in any superlevel. By using the definition of Zi 
and hi = bj this can be expanded to the usual non-LTE partition function 

oo 

Qnon-LTE = £ E^e"^ + £ b i9i ^^ = £ b l9i ^^ (20) 
I t€l igl i=0 

With these relations it is possible to obtain LTE and non-LTE occupation numbers from the 
departure coefficients. If we had used the definition of departure coefficients introduced by Wijbenga 
& Zwaan (1972) the LTE occupation numbers could be obtained simply by solving the Saha- 
Boltzmann equations for the complete molecule. However, we use the definition by Menzel & Cillie 
(1937) and, therefore, it is necessary to use the partition function (19). 

Generally, the partition function is an infinite sum as seen in equation (20) which mathemati- 
cally does not converge. Physically, this sum is "cut off" at some level which can never be populated 
due to the presence of other neighboring particles (see e.g. Mihalas 1978). If the number of levels 
for which data are available is large enough, it is sufficient to set <5 CO rr = and circumvent the 
problem of deciding how to cut off the infinite LTE partition function. 



3.3. The radiative rates within the superlevel formalism 

The emission rate between two superlevels is the sum of all emissions n u {n1 /n*)i? u ; between 
all relevant actual levels. Inserting (n^n^nj) / (nun^n^) and expressing R u i by the cross section 
a e ^ of individual lines (equation (2)) yields : 

n* nfjhcJo f^ nu n* L n* u \ A 5 J 

leL ueu 

and when using equations (13) and (16) 

ts, < n v hc Jo ti Zl v a / 



ueu " " " u lei 

leL ueu 



uu^Rul (21) 



n v 



where 



Rul - ^^°°^i(^ + J A (A))e-^T AdA (22) 

*Wt - zZ a ^f e ~ El/kT ( 23 ) 

leL L 

ueu 

Now, the rate coefficient (22) has the same form and properties (see equation (21)) as the rate 
coefficient (3) for "normal" levels. 
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Similarly, the absorption rate between two superlevels is the sum over all relevant absorptions 
niRiu and becomes 

( A POO ~\ 

n L R LU = n L l—J^ af^j.J\(\)\d\> (24) 

where 

abs _ abs 9l —E,/kT / ri 

tt LU = 2^ a lu z~ e ( 25 ) 

ueu 

Again, the definition of Rlu has the same form as the equivalent expressions (2) for actual levels. 

Note that equations (23) and (25) only differ in the choice of a cm or a abs and that in the 
case of complete redistribution a'fjl and af^ are identical also for superlevels. That means that a 
computer program needs to keep track of only one quantity which is less memory or time consuming 
than keeping track of the sum over the absorption cross sections and the sum over the emission 
cross sections (note that the cross sections need to be known for every depth point of an atmosphere 
and for every transition). 

With the two radiative rates (22) and (24) it is possible to calculate all the radiative rates 
needed to solve the rate equation (4). The rates have to be obtained via wavelength integration 
during the wavelength loop (as described in section (2)) and depend on all actual transitions between 
two different superlevels. Since all actual transitions have to be calculated this does not seem to 
simplify the original problem of dealing with millions of lines. However, phoenix is already designed 
to calculate as many transitions as desired very efficiently, as explained in section (2) by means of 
dOS. Therefore, it is straight forward to calculate the radiative rates "as exactly as desired". The 
computer program simply has to assign each actual transition to the super transition it belongs 
to. Note also, that the purpose of the superlevels in this approach is to reduce the number of 
levels in the rate equations not the number of lines calculated. This keeps the accuracy of the rate 
coefficients and of the absorption coefficients as high as possible. Only the population densities will 
be approximated. It also keeps the resulting spectrum as accurate as possible. 

The rate equations are solved by operator splitting techniques as described in section (2). Note 
that the operator technique only relies on the \E' operator and the E operator. The ^ operator is 
derived from the A operator and the E operator is derived from the emissivity. The A operator 
obviously does not depend on the superlevel formalism. The E operator depends on the opacity, 
but as already explained, the treatment of the opacity is not changed by the superlevel formalism. 

3.4. The absorption and emission coefficients 

As mentioned above, the only purpose of the superlevels is to keep the system of the rate 
equations small. There is no need to keep the number of transitions small. Therefore, every 
absorption or emission coefficient is calculated between actual levels. As explained above, phoenix 
uses lists of spectral lines as input data. The gf values (and other data) for the systems treated 
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with superlevels can be (and are) included in the same input data. Therefore, the cross sections ai u 
are available for every individual transition at every wavelength point selected by dOS as explained 
above. Care has to be taken when entering the number density into the absorption coefficient: 
equation (13) has to be applied to obtain the number density of a particular actual level out of the 
number density of the superlevels. 

That means that, as explained for the radiative rates, our computer code calculates "ex- 
act" cross sections a^u for superlevel transitions with complicated "superline profiles" and does 
not require opacity distribution functions for the supertransitions. The advantage regarding the 
absorption and emission coefficients and the spectrum synthesis is that we can calculate a very 
accurate spectrum at any desired resolution within our superlevel formalism. 



3.5. Collisional and continuum rates within the superlevel formalism 



The collisional excitation rate from one superlevel to another superlevel is the sum of all 
excitation rates between the relevant actual levels and becomes 



n L C LU = n L {Y J f L e- El/kT Ci\. 



ieL 

u£U 



(26) 



Similarly, the rates for collisional de-excitation between two superlevels are : 



nuC UL = njj < 



YJL e -E l/k T c 



ieL 



(27) 



With the definitions (26) and (27) and the relation (16), it follows immediately that 



Cul Zl 



n. 



(28) 



which is the equivalent to the Boltzmann relation for the collisional rate coefficients for actual 
levels. 

With equations (26), (27) and (28) it is now possible to calculate the collisional rates n/C/j 
needed for the rate equation (4) from known collisional cross sections and coefficients Ciu- 

Ionization and recombination processes can be included very similarly. For photoionization 

'4tt ro ° 



" lBl ' "'^T7-J u n '''- /AiA;A,1A 



where 



„ ion 

a L 



ion 9l -Ei/kT 

Z L 



(29) 
(30) 
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and for radiative recombination 



where 

a recomb = J^rccomb ^_ e ~E t /kT ^ 



lei Zl 



For the collisional ionization 



^c = n L |g|-e-^Q c | (33) 

and for the collisional recombination 

n c C ci = n c {|£|e-^Q c } (34) 



However, usually the most important continuum processes of molecules are dissociation and 
molecular recombination. Including these processes requires accounting for all particles that are 
involved in a molecular production and destruction chain, i.e. all constituent atoms and all possible 
molecules that can be constructed out of these atoms via all possible chemical pathways. Such a 
NLCE (non local chemical equilibrium) treatment is beyond the scope of this work since it requires a 
formalism beyond the superlevel formalism. Important molecules with low dissociation energies like 
TiO could be affected by NLCE. However, it can be neglected for molecules that have dissociation 
energies less than those that occur in M dwarf atmospheres, e.g. for CO. 



4. Implementation of the superlevel formalism for CO 

We test our method on carbon monoxide (CO). This diatomic molecule has data available that 
are detailed enough and of good quality (see e.g. Schweitzer 1995). Also, it does not have very many 
lines and levels; and the lines have very distinct bands in the infrared (most prominently at 2.3 /im 
and at 4.3-6 fim). CO not only has reliable data, but is also important in other astronomical fields 
like the earth's atmosphere and the solar atmosphere. 

Other molecules like TiO or H2O have either data which are not detailed enough, i.e. the level 
data cannot be extracted, or the data are most likely incomplete and not very accurate (Allard et al. 
1994; Allard et al. 1997). Furthermore, CO is a relatively small system. With only a few thousand 
levels (about a factor of 1000 less than TiO or H2O) it is only slightly larger than the largest atomic 
system currently calculated in non-LTE. Therefore, we will compare the superlevel calculations with 
a direct non-LTE calculation. Thus, we can quantify the accuracy of the superlevel treatment by 
direct comparison to an "exact" calculation. 
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4.1. Available Data 

4-1.1. Molecular line and level data 

The necessary level data, i.e. the excitation energies, and the quantum numbers as well as the 
necessary line data, i.e. the transition wavelengths, the transition strengths and the involved levels, 
were extracted from the line list by Goorvitch et al. (Goorvitch 1994; Goorvitch & Chackerian, Jr. 
1994a, b). The line list is detailed enough to uniquely assign the upper level to each transition. 

For simplicity, only 12 C 16 was considered. It contains 3623 levels and 19203 transitions for the 
electronic ground state. The first excited electronic state has a very high energy and all transitions 
to this state are in the UV. Therefore, these transitions are not included in the available line list 
as these states and transitions are unimportant in M dwarfs. 



4-1.2. Collisional cross sections 

For collisions with electrons the formula given in Allen (1973) has been used for simplicity. 
This is a very rough approximation but electron densities are very low in M dwarfs and as already 
noted by Thompson (1973) or Ayres & Wiedemann (1989), e.g., electron collisions with CO are 
negligible even for higher electron temperatures than in M dwarfs. 

For collisions with atomic and molecular hydrogen and with helium we used the cross sections 
collected by Ayres & Wiedemann (1989) and expressed in terms of 

O , 0^0-19 ^(^-0-069^/3) 

n ul - 4.2 x 10 0(1 - e -0) ( } 

where (3 = E u [/kT and Q u i = R u i/n x . The subscript x stands for the collisional partner. A x and 
B x are also taken from Ayres & Wiedemann (1989). Explicitly, A x is 3, 64 and 87 for H, H2 and 
He, respectively, and B x is 18.1, 19.1 and 19.1 for H, H2 and He, respectively. These are values 
originating from the results of Glass & Kironde (1982) for H and H2 and the results of Milikan 
(1964) for He. 



4-1.3. Data for continuum processes 

CO has a dissociation energy of 11.09 eV and an ionization energy of 14.01 eV. These energies 
are much too high to be achieved in M dwarf atmospheres (10 eV correspond to 77 200 K or 1240 A). 
Also, no data for detailed cross sections for ionization or dissociation are available. Therefore, the 
molecule is treated in "quasi-internal non-LTE", i.e., our calculations can only determine the 
population densities of the energy levels of CO, not the number density of CO. That means we 
neglect NLCE and the dissociation and ionization transitions are accounted for, but only with 
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negligible transition rates and dissociation processes are not accounted for in the atomic equation 
of state. This will be improved in future work. 

4.2. The division into superlevels 

Energy level diagrams of Carbon monoxide are shown in Figs. (1) and (2). The energy 
levels have been grouped into vibrational quantum numbers and the transitions have been omitted 
for clarity. When dividing the levels into superlevels we have to assure that the coupling within 
one superlevel is very strong. Then we can be certain that the requirements for the superlevel 
approximation are met and that we can assume that we have a thermal population within one 
superlevel. The following possible divisions into superlevels have been considered : 

4.2.1. Model A 

An obvious splitting into superlevels is the grouping by constant vibrational quantum number, 
assuming strong coupling between the rotational states for one vibrational quantum number. This 
assumption is very good since the purely rotational transitions have very long lifetimes. The 
Einstein A values for purely rotational transitions are of the order year" 1 . Therefore, collisions will 
occur much more often and rotational thermalization can be assumed. 

We will call this model "Model A". It corresponds to vibrational LTE as used by Kutepov 
et al. (1991). In Fig. (1), the vertical lines mark the boundaries for Model A. 

However, this grouping results in a large overlap in energy between all superlevels and a wide 
range of energies affects one superlevel. This may make the thermalization assumption invalid. 

4.2.2. Model B 

The second superlevel model selects superlevels only by their energy. The energy boundaries 
are defined by the rotational ground states for each vibrational state. The original authors of 
the superlevel idea (Anderson 1989; Dreizler & Werner 1993; Hubeny & Lanz 1995) also used the 
energy as a selection criteria for their atomic systems. At face value, the pure energy criteria does 
not seem to meet the requirements for the superlevel treatment to be valid. There is, a priori, no 
strong coupling for levels with similar energies. For atoms, this criterion has been a matter of pure 
feasibility. In our case, however, each superlevel defined by the energies of the actual levels will 
still be dominated by the levels with one (the highest) vibrational quantum number. As explained 
above, it is a good approximation to assume strong coupling for constant vibrational quantum 
number. If all the other levels included in such a superlevel that do not have the same vibrational 
quantum number are not important for balancing the rate equations, they do not have to couple 
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strongly to the dominating levels. 

This model we will call "Model B" . Such a superlevel definition can be seen in Fig. (2) indicated 
by the horizontal lines which are the superlevel energy boundaries. 

The vibrational level with the highest quantum number and the other levels with corresponding 
energies would produce a huge and wide spread superlevel. This superlevel was divided in such a 
way that the resulting superlevels contain roughly the same number of levels per superlevel as does 
the average superlevel. 

4.2.3. Model C 

The third superlevel model combines both criteria from Model A and Model B. That means 
we select all levels by their vibrational quantum number and by their energy. That way we have 
the advantages of both methods. The internal coupling can be assured by purely rotational tran- 
sitions which only can occur via collisions (as explained for Model A) and, in addition, there is 
no energetical overlap between the superlevels. Furthermore, this model represents a step between 
Models A or B and a direct solution without superlevels. This model has only a factor of ten less 
superlevels than actual levels. Our model contains 350 superlevels for this method compared to 24 
and 27 for Model A and Model B, respectively. 

We call this model "Model C". It can be visualized by plotting Fig. (2) over Fig. (1). The 
resulting rectangles represent the superlevel boundaries. 

4.2.4. Model Z 

The last model we consider is direct non-LTE using all available levels and lines. The 3623 
levels and 19203 lines of the CO model are within the limits of a modern workstation. This model we 
call "Model Z" . Note, that the superlevel formalism reduces to "regular" non-LTE when the sums 
are carried out over only one term. Therefore, regular non-LTE can easily be treated within our 
superlevel formalism. We will use the results from this model to evaluate our superlevel formalism. 

4.3. Numerical implementation 

The non-LTE implementation for molecules is implemented in a manner very similarly to 
the already existing non-LTE implementation for atomic species as explained in section (2). The 
necessary level data for all actual individual levels are obtained from the input line list. The 
absorption coefficients for individual CO lines are calculated with the non-LTE occupation numbers 
via equation (13). For every wavelength point the cross sections between superlevels are obtained 
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as explained in section (3.3) and are integrated over the spectrum. 

In order to select a certain superlevel model we need only change the input files. The computer 
code is independent of assumptions about the model molecule. This makes it very easy to change 
between the various models. 

The rate equations are solved as described in section (2) using the operator splitting technique. 
Convergence of the departure coefficients b, L is obtained simultaneously with the iteration of the 
temperature structure. 

4-3.1. The accuracy of the rates 

The convergence properties of the rate operator method depends strongly on the accuracy of 
the rates. The accuracy of the rates is determined by the accuracy of the integration of the radiative 
rates. This has been achieved by a carefully selected wavelength grid and by correcting the first 
order error via normalization of the line profiles used for the radiative rates. 

The wavelength grid on which the spectrum is calculated has to include all transitions between 
all superlevels. For the accuracy of the rates it does not have to include the transitions that occur 
within one superlevel. The wavelength grid has been optimized to sample the molecular line profiles 
with as few wavelength points as necessary. Otherwise, the number of necessary wavelength points 
would be at least a multiple of the number of individual lines. For molecules with millions of lines 
that would increase the number of wavelength points to be calculated (which increases the CPU 
time linearly) to an unnecessarily large number. 

4.3.2. LTE tests 

In order to test the stability of the code, several tests have been performed. All these tests 
take advantage of the fact that LTE has to be restored in limiting cases. LTE is restored when the 
atmosphere is dominated by collisions or if the radiation field is described by the Planck function. 
LTE becomes exact when there are only collisional rates and the radiative rates are all zero (Cjj 7^ 0, 
Rij = 0) or when the source function is exactly the Planck function (S v = B v ). Both conditions 
can easily be applied artificially in the computer code. For the second test condition to be very 
clean and accurate, the collisional rates are set to zero (i.e. S v = B u , Cjj = 0). 

All tests have been performed with a model for T e fj=2700, log(g)=5.0 and solar metallicity. 
The choice of superlevels was "Model B" (i.e. rotational ground states as energy boundaries). The 
choice of the superlevel model does not affect this test because a change of the superlevel model is 
only a change of input data, not a change of the computer code. 

The first test case was a LTE model. The departure coefficients 6$ remained 1.0 to an accuracy 
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of about 10 for several iterations for both test conditions. 

The second test case was a LTE model, but the departure coefficients were artificially changed 
under the restriction of particle conservation. That means some fraction (10-50 state were put into 
another excited state. This yields b^s of 0.5 or 0.9 for the de-populated state and frj's of up to ten 
or more for the over-populated state. The LTE situation was restored within only one iteration 
with an accuracy of about 10 for both test conditions as well. 

4-3.3. Numerical stability 

The collisional rates couple all levels with each other. The energy differences (and therefore, 
the rates) between two levels can vary by several orders of magnitudes when comparing arbitrary 
pairs of levels. That has important numerical complications for Model A. In Model A, the individual 
transitions that make up one particular supertransition do have largely varying energy differences 
and rates. Since a collisional rate between two superlevels is the sum of the many individual rates 
we had to assure that the numerically correct sum got evaluated. Otherwise, the collisional rates 
did not obey the Boltzmann relation (28) anymore and the LTE tests described above would fail. 

This is a disadvantage of Model A. All other models do not have this problem since the energy 
differences of the transitions that make up the supertransition are of the same order of magnitude. 
In particular, Model C can be regarded as Model A which has been modified in order to decrease 
the energy differences which then results in numerical stability. For future molecules we will have 
to assure numerical stable summation if possible or otherwise use a different superlevel model. 

4.4. Results for a test molecule 

In the following we present results of a test molecule strongly affected by non-LTE effects. 
We used the CO molecule without accounting for collisions with H, H2 and He. This is a very 
unrealistic model. However, it will be in non-LTE in cool atmospheres. Therefore, we used it to 
test our method in the non-LTE case. 

The models in this section are fully converged and have T c ff=2700 K, log(g)=5.0 and solar 
metallicity. 

4- 4-1- Behavior of the departure coefficients 

The departure coefficients for a converged model can be found in Fig. (3) for Model A, in 
Fig. (4) for Model B, in Fig. (5) for Model C and in Fig. (6) for Model Z. As can be seen, the bi 
structures are very similar for all models. The onset of non-LTE effects start in about the same 
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depth of about r = 10 2 for all superlevel models. 

Generally, it is not expected that the bi distribution is the same for all superlevel models since 
the superlevel models differ. However, Models A and B are expected to be dominated by the same 
levels. The similarity of Figs. (3) and (4) can very well be explained when we assume that Models A 
and B are dominated by the same actual levels. This is also supported by the results from Model C 
which can be seen in Fig. (5). The bright lines in Fig. (5) correspond to the superlevels that have 
the lowest energies for a given vibrational quantum number, i.e. contain the dominating levels for 
Models A and B. As one can see, the bright bi structure in Fig. (5) is very similar to Figs. (3) and 
(4). The strongest argument, however, comes from our direct non-LTE calculation. As can be seen 
in Fig. (6), the levels with the lowest rotational quantum numbers have a very similar bi structure 
as shown in Figs. (3) and (4). That means that, indeed, the superlevels are dominated by the 
levels with the lowest rotational quantum numbers (or equivalently the lowest energies) for a given 
vibrational quantum number and all our superlevel models produce consistent bi structures. 



4- 4-%- Effects on the spectrum 

Detailed high resolution spectra have been calculated in the spectral regions where CO is most 
prominent, i.e. between 2 and 2.4 ^um and between 4.3 and 6 /xm. The Av = 2 band between 2 
and 2.4 /im did not show any significant changes, whereas the Av = 1 band between 4.3 and 6 /ira 
did. 

A high resolution spectrum of a region in the Av = 1 is shown in Fig. (7) for an LTE model and 
non-LTE models with the various superlevel models discussed here. As can be seen, all the non- 
LTE models, including the direct non-LTE calculation, produce deeper lines and are all identical 
within the accuracy of the calculations and the plots. Fig. (7) is only an example of that band. We 
note that all CO lines are deeper for our non-LTE calculations and, therefore, the whole Av = 1 
band is deeper. 

It is important that the spectra are the same for every superlevel model that is chosen under 
the correct physical assumptions. If we want to use the superlevel approximation we have to make 
sure that the physical results are independent of the approximation. But most importantly, the 
direct non-LTE and our superlevel models produce the same result for the spectrum. 

The presence of non-LTE effects in the Av = 1 band also means that this band forms in 
regions where deviations from LTE are present, i.e. outward of r ~ 1CP 2 . The Av = 2 band on the 
other hand has to form in regions deeper than r ~ 10 -2 . This explains that no non-LTE effects 
are visible in the spectral region between 2 and 2.4 fim. 
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4.5. Results for CO 

All models presented here have T e fj=2700 K, log(g)=5.0 and solar metallicity. This roughly 
describes a typical M 8 dwarf. The very low effective temperature has been deliberately chosen 
since non-LTE effects are expected to be largest for low temperatures. The models have been 
calculated using the superlevel models described above and are fully converged with regard to the 
departure coefficients and to the temperature structure. 

4-5.1. Behavior of the departure coefficients 

The departure coefficients for a realistic CO calculation are shown in Figs. (8) (9) (10) and 
(11) for Models A, B, C and Z, respectively. The plots show that CO stays in LTE when accounting 
for collisions with H, He and H2. The levels that deviate significantly from LTE correspond to the 
levels with very high excitation energies and have negligible effects. 

It is very important to note, that we can confirm the previous calculations that found LTE for 
CO not only for approximative calculations but even if we do a direct non-LTE calculation. 

4-5.2. The rate coefficients 

We also compared the collisional rate coefficients due to collisions with electrons, H, He and H2 
with the radiative rate coefficients. As a typical example we show the de-excitation rate coefficients 
-R75 for Model B in figure (12). These are the de-excitations between the superlevels defined by the 
vibrational quantum numbers 7 and 5 and the respective energy bands. As can be seen the rate 
coefficients due to electronic collisions are the smallest collisional rate coefficients. The radiative 
rate coefficients are generally smaller than the collisional rate coefficients, which is pushing the 
system in LTE. 

We also compared the density independent parts of the rate coefficients Q, for the different 
collision partners. As can be seen, the electronic fi's are the strongest, followed by H. However, the 
dominating factor for the resulting rate coefficients are the respective densities. The two largest 
rate coefficients are the one for H because of its large cross section and H2 because of its large 
number density. 

The cross sections for collisions of H, He and H2 with CO are known to be uncertain (Ayres & 
Wiedemann 1989). However, our results indicate that the rate coefficients for H, He and H2 would 
have to be several orders of magnitudes smaller to see non-LTE effects for CO. This is in good 
agreement with all the previous studies. 
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5. Conclusions 

We have developed a method to treat molecules in non-LTE. This method uses the superlevel 
formalism which reduces the number of levels for a model molecule by a factor of hundred or more 
and thereby reducing the rank of the system of rate equations. The important difference with 
respect to other implementations of superlevels — like in the work of Anderson (1989), Dreizler 
&; Werner (1993) or Hubeny & Lanz (1995) — is the treatment of line opacities. The original 
accuracy of opacity treatment implemented in the atmosphere code phoenix is still preserved by 
taking advantage of the direct opacity sampling (dOS). This treats the line opacity as exactly as 
possible and allows us to also use this also for the radiative line rates. 

We used CO as an example and showed how it is possible to create superlevel models that 
approximate the direct calculations very well. The most important result we can draw from CO is 
the fact that we successfully divided the CO molecule into a superlevel model molecule. The results 
we obtained are independent of the choice of the superlevel model and are identical to the results 
a direct non-LTE calculation of CO which is the largest non-LTE calculation of a single species 
ever. Therefore, we demonstrated that our superlevel method is a powerful and accurate method 
to calculate molecular non-LTE problems. For physical reasons Model A has to be recommended 
as long as the numerical stability is guaranteed. 

For CO itself we can confirm the results from the previous studies by, e.g., Ayres & Wiedemann 
(1989, and references therein) that CO is in LTE in cool stellar atmospheres. 

In future work we will apply our method to the most important, but much more complicated 
molecules TiO and H2O. The results from CO suggest that it does not matter much which superlevel 
we choose as long as it is physically reasonable and numerically stable. Therefore, it might be 
tempting to use Model A for TiO and H2O since Model A is the most physical and gives correct 
results. However, we will investigate several superlevel models also for TiO and H2O which will be 
similar to Model B and C or variations thereof. Thus we can assure the physical correctness, the 
numerical stability and the necessary decrease of the system of rate equations. 

For molecules other than CO it will be important to also include good cross sections for 
continuum processes. Continuum processes for CO could be ignored, since CO is the most stable 
molecule. However, when calculating molecules like TiO or H2O it will be important to allow for 
ionization and recombination, and, most importantly, for dissociation and molecule generation. In 
many cases dissociation will dominate over ionization. This will couple the atomic equation of 
state and the molecular equation of state by sources and sinks of atoms and ions and production 
chains for molecules have to be calculated simultaneously. Therefore, the logical consequence is 
to develop a formalism for non local chemical equilibrium (NLCE). NLCE will be important for 
molecules with low dissociation energies like TiO. Non-LTE effects on Ti I lines have already been 
presented in Hauschildt et al. (1997a), therefore, it will be necessary to treat TiO, Ti I and possibly 
O I simultaneously in non-LTE. 
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Whereas pure ionization can be done within the framework of the superlevel formalism, the 
inclusion of NLCE will require further technical developments. NLCE will not affect the rate 
equations by increasing the number of levels but by coupling to other species. When developing 
a NLCE formalism, it will be based on the superlevel formalism, the superlevel formalism for 
molecules as presented in this work for very stable molecules like CO. 

Finally, since it is now possible to calculate CO with a lot of levels in non-LTE, these methods 
can now be applied to fields other than M dwarf atmospheres. One application will be to address 
the problem of the strong CO bands in the solar spectrum as already mentioned. 
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Fig. 1. — Energy level diagram for the electronic ground state of CO. The bound-bound transitions 
have been omitted for clarity. The vertical lines mark the superlevel boundaries for Model A (see 
text for details). 
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Fig. 2. — Energy level diagram for the electronic ground state of CO. The bound-bound transitions 
have been omitted for clarity. The horizontal lines mark the superlevel boundaries for Model B 
(see text for details). 
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Fig. 3. — The departure coefficients for a converged modei with temperature correction. Selection 
method is "Model A" , i.e. the vibrational quantum number. The temperature scale on the top is 
the electron temperature at the respective depth. 



-27- 



Temperature in K 
1583 1639 1826 2147 2848 3426 




10" 10" 10" 10 w 10 

Optical depth 



Fig. 4. — The departure coefficients for a converged model with temperature correction. Selection 
method is "Model B", i.e. the energy of the rotational ground states. The temperature scale on 
the top is the electron temperature at the respective depth. 
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Fig. 5. — The departure coefficients for a converged model with temperature correction. Selection 
method is "Model C" , i.e. the energy and the vibrational quantum number. The temperature scale 
on the top is the electron temperature at the respective depth. The bright lines correspond to the 
superlevels that have the lowest energies for a given vibrational quantum number. 
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Fig. 6. — The departure coefficients for a converged model with temperature correction. The model 
is "Model Z", i.e. direct non-LTE. The temperature scale on the top is the electron temperature 
at the respective depth. The levels with the lowest rotational quantum numbers are indicated in 
the figure. 
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Fig. 7. — High resolution model spectra for an arbitrary region in the Av = 1 band between 4.3 and 
6 /um. The LTE spectrum and the non-LTE spectra with the different models for the superlevel are 
indicated in the figure. All non-LTE spectra lie on top of each other. This spectrum is calculated 
with a very high resolution to demonstrate the effects on the individual lines. 
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Fig. 8. — The departure coefficients for a converged model with temperature correction. Selection 
method is "Model A" , i.e. the vibrational quantum number. The temperature scale on the top is 
the electron temperature at the respective depth. 
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Fig. 9. — The departure coefficients for a converged model with temperature correction. Selection 
method is "Model B", i.e. the energy of the rotational ground states. The temperature scale on 
the top is the electron temperature at the respective depth. 
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Fig. 10. — The departure coefficients for a converged model with temperature correction. Selection 
method is "Model C" , i.e. the energy and the vibrational quantum number. The temperature scale 
on the top is the electron temperature at the respective depth. The bright lines correspond to the 
superlevels that have the lowest energies for a given vibrational quantum number. 




Fig. 11. — The departure coefficients for a converged model with temperature correction. The model 
is "Model Z", i.e. direct non-LTE. The temperature scale on the top is the electron temperature 
at the respective depth. The levels with the lowest rotational quantum numbers are indicated in 
the figure. 
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Fig. 12. — The collisional de-excitation rate coefficients, the radiative emission coefficients and the 
Q factors for the transition between superlevels number 7 and 5 in Model B (superlevel number 
is the superlevel with the lowest excitation energies). The collisional partners are indicated in the 
figure. The top part of the figure are collisional rate coefficients, the bottom part are the f2 factors. 
The units are phoenix internal arbitrary units, but comparable among each other. The O factors 
have to be multiplied with the densities of the collisional partners to obtain the rate coefficients. 



